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LONG-TERM  GOALS 

The  long  term  goal  of  this  project  is  to  develop  inversion  approaches  that  enable  the  estimation  of 
refractivity  profdes  and  the  associated  uncertainty.  Furthermore,  we  will  develop  methods  for 
mapping  the  refractivity  parameters  and  their  associated  uncertainty  into  propagation. 

OBJECTIVES 

The  objective  of  this  proposal  is  the  development  of  parametric  approaches  for  the  inversion  of  radar 
clutter  data  to  estimate  atmospheric  refractivity  over  land  and  sea.  Refractivity  inversion  algorithms 
using  land  clutter  with  and  without  radar  phase  will  be  developed.  Over  the  sea  we  will  further  develop 
our  approaches  to  carry  out  estimation  in  time  and  space  of  the  refractivity. 

APPROACH 

We  have  considerable  experience  in  carrying  out  refractivity  estimation  from  ocean  clutter  data 
[Gerstoft  et  ah,  2003a,  2003b,  Gerstoft  et  ah,  2004;  Rogers  et  ah,  2004].  Little  has  been  done  to 
indicate  the  quality  of  the  solution  for  each  parameter,  either  with  the  variance  of  the  parameter 
estimate  or  preferably  the  complete  a  posteriori  distribution.  We  have  already  done  much  work  on  this 
in  an  ocean  acoustic  context,  but  this  has  not  been  explored  in  our  refractivity  from  clutter  (RFC) 
processing  to  date.  This  will  entail  developing  likelihood  formulations  and  importance  sampling 
algorithms.  This  inversion  approach  will  show  the  information  content  in  the  data,  the  importance  of 
each  parameter,  and  the  quality  of  the  inversions.  Another  related  topic  is  that  in  RFC  inversions,  we 
commonly  invert  each  data  block  independently.  When  these  inversions  are  close  in  time  (i.e., 
successive  looks  in  time  at  the  same  azimuth)  or  space  (i.e.  adjacent  azimuths),  it  should  be  beneficial 
to  use  the  results  of  the  previous  inversion  as  a  starting  condition  for  the  next  inversion.  A  natural 
framework  for  this  is  a  Bayesian  approach  where  the  posterior  of  the  last  inversion  becomes  the  prior 
for  the  current  inversion.  For  this  investigation,  the  Bayesian  approach  will  be  implemented  using  a 
Metropolis-Hastings  Gibbs  sampler. 

WORK  COMPLETED 

A  method  for  estimation  of  the  radio  refractivity  from  radar  clutter  using  Markov  Chain  Monte  Carlo 
(MCMC)  samplers  with  a  likelihood-based  Bayesian  inversion  formulation  has  been  introduced.  This 
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approach  enables  us  to  obtain  full  ^-dimensional  posterior  probability  distributions  for  the  unknown 
parameters  as  well  as  the  maximum  likelihood  solution  itself  [Yardim  05a,  05b,  06]. 
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Figure  1  Clutter  map  from  Space  Range  Radar  (SPANDAR)  at  Wallops  Island^  VA. 
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Figure  2  Tri— linear  M-profile  and  its  corresponding  coverage  diagram. 
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Figure  3  The  4-parameter  tri-linear  M-profile  model  used  in  this  work. 


RESULTS 

An  accurate  knowledge  of  radio  refractivity  is  essential  in  many  radar  and  propagation  applications. 
Especially  at  low  altitudes,  radio  refractivity  can  vary  considerably  with  both  height  and  range,  heavily 
affecting  the  propagation  characteristics.  One  important  example  is  the  formation  of  an 
electromagnetic  duct.  A  signal  sent  from  a  surface  or  low  altitude  source,  such  as  a  ship  or  low-flying 
object,  can  be  totally  trapped  in  the  duct.  This  will  result  in  multiple  reflections  from  the  surface  and 
they  will  appear  as  clutter  rings  in  the  radar  PPI  screen  (Fig.  1).  In  such  cases,  a  standard  atmospheric 
assumption  with  a  slope  of  modified  refractivity  of  0.1 18  M-units/m  may  not  give  reliable  predictions 
for  a  radar  system  operating  in  such  an  environment. 

Ducting  is  a  phenomenon  that  is  encountered  mostly  in  sea-borne  applications  due  to  the  abrupt 
changes  in  the  vertical  temperature  and  humidity  profiles  just  above  large  water  masses,  which  may 
result  in  an  sharp  decrease  in  the  modified  refractivity  (M-profile)  with  increasing  altitude.  This  will, 
in  turn,  cause  the  electromagnetic  signal  to  bend  downward,  effectively  trapping  the  signal  within  the 
duct.  It  is  frequently  encountered  in  many  regions  of  the  world  such  as  the  Persian  Gulf,  the 
Mediterranean  and  California.  In  many  cases,  a  simple  tri-linear  M-profile  is  used  to  describe  this 
variation.  The  coverage  diagram  of  a  trapped  signal  in  such  an  environment  is  given  in  Fig.  2. 

As  detailed  in  the  paper  [Yardim  at  2006]  we  have  developed  a  likelihood  based  Markov  Chain  Monte 
Carlo  (MCMC)  sampler,  which  we  have  compared  to  a  classical  genetic  algorithm  and  an  exhaustive 
grid  sampling.  The  exhaustive  grid  sampling  is  only  possible  for  search  spaces  of  small  dimensions. 

In  previous  work,  genetic  algorithms  (GA)  and  Markov  chain  Monte  Carlo  (MCMC)  samplers  were 
used  to  calculate  the  atmospheric  refractivity  from  returned  radar  clutter.  Although  GA  is  fast  and 
estimates  the  maximum  a  posteriori  (MAP)  solution  well,  it  poorly  calculates  the  multi-dimensional 


integrals  required  to  obtain  means,  varianees  and  underlying  posterior  probability  distribution 
funetions  (PPD)  of  the  estimated  parameters.  Aecurate  distributions  and  integral  ealculations  can  be 
obtained  using  MCMC  samplers,  such  as  the  Metropolis-Hastings  (M-H)  and  Gibbs  sampling  (GS) 
algorithms.  Their  drawbacks  are  that  they  require  a  large  number  of  samples  relative  to  optimization 
techniques  such  as  GA  and  become  impractical  with  increasing  number  of  unknowns. 

A  hybrid  GA-MCMC  method  based  on  the  nearest  neighborhood  algorithm  (NA)  is  implemented.  It  is 
classified  as  an  improved  GA  method  which  improves  integral  calculation  accuracy  through 
hybridization  with  a  MCMC  sampler.  Since  it  is  mainly  GA,  it  requires  fewer  forward  samples  than  a 
MCMC,  enabling  inversion  of  atmospheric  models  with  a  larger  number  of  unknowns. 

MCMC  samplers  are  much  faster  than  exhaustive  search,  but  still  need  much  more  forward  model  runs 
than  GA.  Therefore,  it  will  be  very  desirable  to  have  a  GA-MCMC  hybrid  method  that  combines  the 
speed  of  GA  with  the  accuracy  of  MCMC.  Due  to  its  inherent  properties,  Gibbs  sampler  (GS)  is  the 
best  MCMC  algorithm  for  such  a  hybrid  method.  Hence,  the  hybrid  method  will  be  called  the  GA-GS 
hybrid.  The  method  consists  of  three  distinct  sections: 

1.  GA:  Run  a  classical  GA,  minimizing  the  misfit,  and  save  all  the  populations  and  the  misfit 
values  of  all  generations. 

2.  Voronoi  Cells  and  Approximate  PPD:  Using  the  GA  samples  and  their  likelihood  values 
construct  Voronoi  cells  (see  next  section)  around  each  GA  point,  assigning  the  likelihood  of 
each  GA  point  to  all  points  inside  its  own  Voronoi  cell.  This  will  result  in  an  approximate  PPD. 

3.  GS:  Run  a  fast  GS  on  the  approximate  PPD  instead  of  the  real  one.  Since  the  conditional 
densities  and  the  likelihood  values  required  to  run  GS  is  known  for  the  approximate  PPD,  no 
forward  model  is  needed. 

This  process  simply  is  a  discretization  of  the  original  PPD.  It  will  convert  the  true  analog  PPD  into  a 
digital  one  through  an  A/D  converter.  The  only  difference  is  that,  this  A/D  converter  is  ^-dimensional, 
and  hence,  discrete  levels  are  tiny  n-dimensional  hypercubes.  These  hypercubes  are  called  Voronoi 
cells.  The  shapes  and  sizes  of  these  cells  are  determined  by  the  GA  sample  set.  There  is  only  one  GA 
sample  in  each  cell  and  there  does  not  exist  any  other  GA  sample,  which  is  closer  to  any  point  inside 
this  hypercube.  Hence,  for  any  boundary  point  between  any  two  adjacent  cells,  the  distances  of  that 
point  to  the  two  closest  GA  samples  are  same.  Due  to  this  feature,  it  is  also  called  the  nearest 
neighborhood  (NN)  method  [Sambridge,  1998].  In  each  cell,  the  likelihood  value  is  assumed  to  be 
constant  with  a  value  of  its  GA  sample.  Therefore,  the  likelihood  of  any  point  anywhere  in  the  entire 
search  space  is  known  and  there  is  no  need  for  any  further  forward  model  runs.  A  very  fast  GS, 
without  any  forward  modeling,  is  used  to  sample  this  approximate  PPD.  The  accuracy  of  the  results 
depends  mostly  on  the  quality  of  the  approximate  PPD,  which  means  that,  GA  should  gather  enough 
samples  from  all  over  the  n-dimensional  search  space  to  allow  the  NN  algorithm  to  construct  a  good 
enough  n-dimensional  mesh,  hence  a  good  enough  approximate  PPD. 

The  method  has  been  carefully  validated  in  Yardim  et  al  [2007].  Here  we  present  inversion  results  of 
clutter  return  from  the  1998  Wallops  experiment  as  also  described  in  Yardim  et  al  [2007].  We 
represent  the  environment  using  a  range  dependent  refractivity  profile  as  given  in  Fig  5.  The 
environmental  posterior  density  is  given  in  Fig.  6(a).  Since  the  full  PPD  is  16-D,  only  1-D  (diagonal 


plots)  and  2-D  (upper  diagonal)  marginal  densities  ealculated  are  given.  Some  of  the  parameters  sueh 
as  m\Q,  m\2  and  m  14  have  a  highly  non-Gaussian  marginals,  while  others  sueh  as  m2,  m3  and  mg  have 
Gaussian- like  features.  The  highly  skewed  1-D  marginals  given  for  m  10  and  w  14  are  eneountered 
frequently  with  the  reffactivity  slope  pdfs.  The  reason  is  that  the  slope  very  rarely  exceeds  values  such 
as  0.3— 0.4  M-units/m  and  usually  is  concentrated  around  values  such  as  0.1 18  M-units/m  (standard 
atmosphere)  and  0. 13  M-units/m.  This  creates  a  sharp  peak  for  the  positive  end  of  the  spectrum  since 
the  negative  slope  values  can  be  in  excess  of  the  -2  M-units/m,  usually  with  a  quickly  decreasing 
probability.  The  result  is  a  pdf  structure  similar  to  the  ones  obtained  here.  In  fact,  Gerstoft  et  al  [2004] 
uses  such  a  pdf  as  prior  density  to  do  importance  sampling. 

The  environmental  statistics  can  be  projected  into  statistics  for  user  parameters  (see  Sec  II  of  Yardim 
et  al  [2007]).  One  typical  parameter  of  interest  to  an  end-user  is  the  propagation  factor  F.  The  results  in 
Fig.  7  are  obtained  from  the  parameter  PPD  in  Fig.  6.  It  shows  the  PPD  for  F  at  ranges  (a)  18,  (b)  40, 
and  (c)  60  km.  Contour  plots  show  the  PPD  of  F  for  height  values  between  0-200  m,  with  the  MAP 
solution  (dashed  white).  Horizontal  lines  represent  the  three  altitudes  analyzed  in  detail  in  the  small 
plots  shown  next  to  the  color  plots.  Comparison  of  plots  at  the  same  range  and  different  altitudes 
reveals  some  important  aspects  of  RFC. 
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Figure  4  Voronoi  cells  for  a  simple  2-parameter  search  space  and  an  approximate 
conditional  1-D  density  for  a  given  cut.  Dots  represent  GA  points  and  squares 
represent  the  GS points  sampling  the  approximate  2-D  PPD. 
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Figure  5  Range-dependent  sixteen  parameter  M-profile  with  four  parameters  per 
20  km.  Vertical  profile  at  any  given  range  is  calculated  by  linear  interpolation  of 
both  the  slopes  and  the  layer  thicknesses. 
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Figure  6:  Marginal  and  conditional  distributions.  (a)l-D  (diagonal)  and  2-D  (upper  diagonal) 
posterior  probability  distributions  in  terms  of  percent  HPD,  for  the  range-dependent  SPANDAR 
data  inversion.  13  parameters  out  of  16  are  shown .  Vertical  lines  in  the  1-D  plots  show  the  GA 
MAP  solution,  (b)  Normalized  error  function  for  various  conditional  planes.  Each  2-D  plot  is 
obtained  by  fixing  the  other  14  parameters  to  their  MAP  values. 


Figure  7  Posterior  probability  densities  for  propagation  factor  F  at  three  ranges:  (a)  18,  (b)  40,  and 
(c)  60  km.  Color  plots  show  the  PPD  of  F  for  height  values  between  0  m  and  200  m  in  terms  of 
percent  HPD,  with  the  MAP  solution  (dashed  white).  Horizontal  lines  represent  the  three  altitudes 
analyzed  in  detail  in  the  small  plots  shown  next  to  the  color  plots  at  heights  180, 100,  and  20  m, 
respectively  from  top  to  bottom.  Vertical  lines  in  the  small  plots  represent  the  values  of  F  at  the 
corresponding  height  and  range  for  the  MAP  solution  (blue  line  with  circles),  helicopter 
measurement  (red),  and  standard  atmospheric  assumption  (black). 


IMPACT/APPLICATIONS 

Knowledge  of  refractivity  profiles  is  important  for  radar  performance  prediction.  Using  the  radar 
clutter  return  to  estimate  refractivity  gives  s  a  real-time  estimate  of  the  refractivity. 

RELATED  PROJECTS 

Refactivity  Data  Fusion  and  Assimilation  (Ted  Rogers,  SPA  WAR):  This  project  is  concerned  with  near 
real-time  techniques  for  inferring  refractivity  parameters  from  radar  sea  clutter. 
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